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We use the conformal approach to numerical relativity to evolve hyperboloidal 
gravitational wave data without any symmetry assumptions. Although our grid 
is finite in space and time, we cover the whole future of the initial data in our 
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calculation, including future null and future timelike infinity. 
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In the articles [p], @, 3 we presented a complete code for solving the conformal Einstein 
equations which will allow us to study many interesting questions about the global struc- 
ture of spacetimes by performing numerical experiments. In the present paper we want to 
demonstrate some of the unique capabilities of this code by a simple example: We calculate 
the time evolution of gravitational wave data without continuous symmetries. The data 
are prescribed on a hyperboloidal, whence spacelike, initial slice extending to future null 
infinity. Starting from this slice, we calculate a conformal spacetime (M, g a i>) by solving 
the conformal time evolution equations. On the region M of the conformal spacetime on 
which the variable Q, the conformal factor, is positive, the metric g a b = Q~ 2 gab defines an 
asymptotically flat solution to Einstein's vacuum field equations. We call (M, g a b) the phys- 
ical spacetime. The boundary of M in M, given by the set {Q = 0}, represents future null 
^ \ infinity, J7" + , and future timelike infinity, i + . 

In our evolution the relevant part of the set {Q = 0} is completely covered by a numerical 
grid. Therefore, embedding the physical spacetime into a larger conformal spacetime implies 
that the determination of the gravitational radiation is a well-defined, gauge ambiguity-free 
procedure and that we avoid any influence of artificial boundaries. It also allows very accu- 
rate determination of the fall-off behaviour of physical quantities near the different infinities. 
Our data are chosen sufficiently close to Minkowski data, so that the solution admits a regu- 
lar point i + in the conformal extension. However, they are also far enough away to produce 
a spacetime which differs significantly from Minkowski space. It has been known theoreti- 
cally for some time that sufficiently weak data should admit a regular point i + at timelike 
infinity || (cf. also Ij) for similar results for weak Cauchy data). However, it is a quite re- 
markable fact that this point can be modelled in a numerical calculation with the precision 
discussed below. 

In the present paper we shall not attempt to discuss the background of the conformal 
approach again. We refer the reader to M and to the recent survey article by J. Frauen- 
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diener JF|. In section II we shall describe the given data. Then we discuss their evolution 
(section |111|) and show, in particular, that we have indeed covered the initial hypersurface 
as well as future null infinity and timelike infinity. 



II. THE INITIAL DATA 

To calculate permissable initial data we use the numerical scheme described in [||], to which 
we refer for details. There, on an initial hypersurface we give a boundary defining function 
0, which will be related to the conformal factor Q and the solution <p of equation ([[]) by 
Q = Q/(j), and a conformal 3-metric h ab . The 3-metric is chosen such that the tracefree part 
of the extrinsic curvature of the surface Q = 0, the initial cut S of null infinity, must vanish. 
Then we solve the Yamabe equation, 

4 n 2 (3) A0 - 4 ft( (3 V a ft) ( (3 V a 0) 

- (3) i?n 2 + 2n (3) An-3( (3 V a fi)( (3 V a fi))0--P0 5 = 0, (1) 

where (3) V, (3) A, and (3) i? are the derivative operator, the Laplace operator, and the Ricci 
scalar associated with h ab , and A; is a positive constant, which we choose to be k = 3. 
The solution of the Yamabe equation, the given free functions Cl and h abl and two more 
free functions, namely the trace k of the conformal extrinsic curvature of the initial slice and 
the Ricci scalar R of the conformal spacetime, define a set of data for the conformal field 
equations. 

Our choices for the free functions are 

n= l -(l-(x 2 + y 2 + z 2 )) (2a) 
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(2b) 



k = (2c) 
R = 0. (2d) 

If is a solution of the Yamabe equation with respect to Cl and h a b, then for any positive 
function 9 the solution of the Yamabe equation with respect to OVl and h a b is 4>V6. Therefore, 
the calculated data depend on the location of S but not on the choice of Cl in the interior, 
and there is no loss of generality, if we choose Cl to be spherically symmetric. 
The 3-metric h ab is obtained by perturbing the xx_ component away from Minkowski data. 
The chosen perturbation has no obvious continuous symmetry, but it is reflection symmetric 
at the planes x — 0, y — 0, and z = 0. For that reason it is sufficient to plot the octant 
{(x, y, z)\x > 0, y > 0, z > 0}, although the calculation has been performed for (x,y,z) 6 
[-1.25, 1.25] x [-1.25, 1.25] x [-1.25, 1.25]. 

The choices of k and R are pure gauge. The function k determines the initial time derivative 
Qq of the conformal factor Q through the relation (7) of 0, the choice of R determines the 
values of certain curvature components (1,1) R a b and E ab and the variable u>, obtained by 
applying the wave operator applied to Q. 

To calculate initial data, we use a spectral bases of 46 x 44 x 44 elements. The maximum 
value of the constraint violation on the initial slice is then of order 10~ 4 . 
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The left graphic in figure [1] shows the metric component h^x on the z = plane for a value 
ofA = l. ~ 




FIG. 1: The initial values for the metric component h xx (left) and the real part of the curvature 
invariant I (right) on the z = plane. Only the physical part of the grid is shown. 

The plane z = is the plane, on which the perturbation of h^x assumes its maximum in 
the physical region, which is approximately 0.025. A priori, it is not clear, whether the 
perturbation is not just a coordinate effect. The right graph in figure H shows the real part 
9ft (/) of the curvature invariant 

/ = $ABCD ¥ BCD 

= n 6 (E ab E ab - B ab B ab + 2 % E ab B ab ) , (3) 

where iPabcd is the Weyl spinor of the physical spacetime, E ab the electric part of the 
rescaled conformal Weyl tensor, B ab its magnetic part, and i = y— 1. Since this curvature 
invariant does not vanish, our data do not correspond to Minkowski space. 
The curvature invariant 9ft (/) assumes its maximum value of approximately 0.00068 at the 
origin. The maximum value for an amplitude A = 2.0 is approximately 0.00270. The 
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curvature invariant dt(I) is therefore approximately proportional to A 2 . There is another 
local maximum at y ~ 0.55 with a value of approximately 0.00061. The physical length of 
the y coordinate line connecting this maximum with the origin is about 1.25. 

III. THE TIME EVOLUTION 

We split the discussion of the time evolution of our data into three parts. The first part 
describes the conformal structure of the spacetime. In the second part we describe how we 
proceed to reconstruct properties of the associated physical spacetime, i. e. the proper time 
of observers. Issues related to gravitational radiation, such as the longtime decay of the 
Bondi mass, are dealt with in the third part. 

To evolve our data, we have to choose five gauge source functions, namely q = \n(N/^/h), 
where N is the lapse and h the determinant of h a b, the three components of the shift iV a , 
and the Ricci scalar R. In our case the simplest choice, 

q = (4a) 
N a = (4b) 
R = (4c) 

is sufficient to cover the entire future of the initial data. 

The size of the next numerical time step is calculated from the Courant-Friedrich-Levy 
condition at each step as described in ||. Corresponding time slices of runs with different 
resolutions do in general not coincide, since the size of the time step depends on the 3-metric 
h a b, which is a variable of our system. 

Simultaneous to the time evolution equation we solve 15 ordinary differential equations 
describing geodesies of the physical metric g a b (cf. Q for numerical details). These world 
lines represent observers in the physical spacetime. Initially the observers are placed at 
coordinate values of (0, 0, 0), (±±, 0, 0), (0, ±§, 0), (0, 0, ±|), and (^L, ^). As initial 
tangent vector we choose the normal of the initial slice and normalise it with respect to the 
physical metric g a b- 

In addition to the observers in physical spacetime we calculate the orbits of 1986 "Bondi 
observers moving along generators of null infinity" . Their orbits define a discretisation of J 
and enable us to calculate radiative quantities (cf. j| for details). The Bondi observers are 
placed at a uniform grid parametrising the initial cut S of J . They are placed at the north 
and the south pole and at 64 x 31 gridpoints covering tp) e]0, 7r[x[0, 27r[. 
An important property of every numerical simulation is an error estimate. The convergence 
analysis of runs with spatial grids of (50 3 ), (100 3 ), and (150 3 ) gridpoints performed on slices 
at t — 0.0, 0.24, 0.49, 0.72, and 0.92 gives an estimate for the absolute maximal error in any 
variable of less than 0.001. This value agrees with the errors obtained when reproducing 
exact solutions with the same code [Q. We have also performed convergence analyses for 
each plot shown. If not explicitely mentioned, the error estimate of the convergence analysis 
suggests, that the error is at most of the order of the line thickness in the plots. 

A. The conformal spacetime 

During the evolution of our data we monitor in particular the behaviour of the conformal 
factor Q. We will have covered the entire physical future of the initial data, if the region 
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{Q > 0} vanishes for one slice, since the physical spacetime is identical with the set {Q > 0}. 
The orbits of the Bondi observers generate the surface {Q = 0}. In figure we plot their 




FIG. 2: The projection of representative null 
generators of J onto the (z, t) plane near the 
focal point. 

orbits near their focal point. 

Due to practical reasons we restrict ourselves to the projection onto the plane {x = 0, y = 0} 
and a representative selection of our 1986 Bondi observers, namely the observers initially 
placed at the 33 points (i?, ip) = (0, 0), (vr/32, 0), . . . , (tt, 0). Plots showing the projections 
to the planes {y = 0, z = 0} and {x = 0, z = 0} and projections of all the other orbits give 
corresponding results. 

It should also be pointed out that the size of a grid cell in the 150 3 run we present is 
At x Ax x Ay x Az pa 0.003 x 0.0 1 7 3 near the focal point. The generators of J clearly meet 
within the volume of one grid cell. The result from a 100 3 run is visibly indistinguishable 
from the 150 3 run, in the 50 3 run the focal point has a slightly smaller t coordinate, namely ~ 
0.996, which is still in excellent agreement, since the size of the 50 3 grid cells is ~ 0.009 x 0.05 3 
near the focal point, which is larger than the deviation of the runs. 

We have continued our calculation beyond the focal point up to t = 1.1 to check regularity of 
the conformal spacetime at the focal point. Since we have completely covered the physical 
future of our initial slice at t — 1 already, it makes no sense to integrate further. Even 
beyond the focal point the conformal spacetime stays regular. Therefore the focal point is 
an excellent candidate for a regular i + . And indeed, as we will see later, physical observer 
will reach this point after infinite proper time. 

Since we have found a complete J + and a regular i + , the constructed spacetime possesses 
qualitatively the same asymptotic structure as the Minkowski space. Quantitatively, there 
are differences in the asymptotics as can be seen from figure |3|, where we plot the real part 
9ft(J) of the conformal curvature invariant 

L — Wabcd W 

= E ab E ah - B ab B ab + 2 i E ab B ab (5) 
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FIG. 3: Time evolution of real part of confor- 
mal curvature invariant / on the positive y axis 
(only the physical region of conformal spacetime 
is shown). 

on the physical portion of the grid. The quantity iPabcd denotes the rescaled Weyl spinor. 
There is a significant amount of rescaled conformal curvature at J + and at i + , indicating 
a non-vanishing fall-off coefficient in the expansion of the corresponding physical curvature 
invariant dt(I) at infinity. 

B. Reconstructing the associated physical spacetime 

When we reconstruct the physical spacetime from the conformal spacetime, we have to con- 
sider two types of quantities. 

The first type consists of those quantities which can be expressed as a regular expression 
in the variables of the conformal field equations. Figure |] shows the time evolution of the 
value of 9ft(i), which is such a quantity, on the y axis. We see that the time evolution of the 
two maxima of the initial data is dominated by a rapid decay towards future null and future 
timelike infinity, although there is another, smaller extremum forming in between the initial 
maxima. This smaller extremum can best be seen in the contour plot on top of the surface 
plot. In the contour plot we have also marked the location of J + . 

With a code performing the numerical integration in physical spacetime we would only be 
able to calculate figure £|, but not figure §. Since the physical curvature invariant decays so 
rapidly to zero (at a conformal time of 0.6 it has already decayed by a factor of order 10000), 
only a very accurate physical code could resolve the fall-off at a large physical time. In the 
conformal picture the decay has been factored out by choosing the rescaled conformal Weyl 
tensor as variable. We calculate the conformal invariant and the conformal factor Q, 
which do not change dramatically during the whole time evolution, and then get the physical 
invariant by multiplying the conformal curvature invariant with the appropriate powers of 
the conformal factor (equations [| and ||). 

A convergence analysis shows that the numerical prediction of the behaviour of 9ft(2) for a 
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FIG. 4: Time evolution of the real part of the 
physical curvature invariant I on the positive y 
axis. 

150 3 calculation is indistinguishable from the 100 3 calculation for the whole calculation up 
to i + . At the last slice of the 100 3 run before i + the curvature invariant has decayed by a 
factor of more than 10 20 . Due to the effect of rounding errors, a physical code working with 
8 byte reals could not resolve the decay over such a large range, regardless of the grid size 
required to achieve the required accuracy. 

The physical metric g a b = VL~ 2 g ah blows up at infinity, since the conformal metric g a b is 
regular everywhere. The second type of quantity consists of those quantities which describe 
physical distances. They, of course, must blow up when approaching i + or J + . 
A typical example is the proper time of observers. We will see in the following paragraphs, 
that the question "For how long can we numerically calculate the measurements of an ob- 
server?" is closely tied to the numerical question "How well can we resolve the neighbourhood 
of the surface {Q = 0}?". 

Figure |^ shows the world line of a representative observer, who is initially placed at 
(x,y,z) = (7^75; 7^75; 27/3-'' * n an ( r = V x<2 + V 2 + z 2 , t) plot. The world line runs into 
i + as expected. 

We can plot the x, y, and z coordinates of the world line of the observer as a function of 
conformal time. If our spacetime were a representation of Minkowski space, and if we were 
to choose the same gauge source functions, the differences x — y, x — z, and y — z would 
vanish for all times, due to symmetry reasons. Figure ^ shows, that this is not the case in the 
computed spacetime. The observer oscillates around the Minkowski orbit — he is pushed 
around by the gravitational wave. The amplitude of the oscillation is small compared to the 
coordinate values. By plotting the difference we make the oscillation visible. 
A plot of the conformal factor along the world line of our observer gives the left graph of 
figure 0. The right graph of that figure shows the vicinity of i + for the 50 3 (x), the 100 3 
(+), and the 150 3 (*) runs. Obviously, the finer the resolution the closer we get to the 
(quadratic) zero of Q. 

When calculating the physical geodesies we also calculate the proper time r as a function 
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Wordline of Observer and Scri 




0.0 0.2 0.4 0.6 0.8 1.0 

Coordinate distance from origin 



FIG. 5: The thin line represents the world line of 
an observer (timelike physical geodesic) starting 
at (x,y,z) = (^=,^=,^=) and going to i+. 
The thick line represents J + . The particular 
behaviour of the world line near i + is due to the 
particular choice of its initial velocity. 




FIG. 6: Oscillation of an observer in coordi- 
nate space: The differences x — y (solid), x — z 
(dashed), and y—z (dotted) as a function of con- 
formal time t which would all be for Minkowski 
data (A = 0). 
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for observer 
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Q for observer near i + 
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FIG. 7: Left: Conformal factor f2 along an observer world line for the entire time evolution. 
Right: The same near i + (right) for a 50 3 (x), an 100 3 (+), and an 150 3 (*) run. 

of the conformal time t. In Figure § we show the result for our observer near i + . In the 50 3 
run we have covered a proper time interval of [0, 300], in the 100 3 run an interval of [0, 350], 
and in the 150 3 run an interval of [0, 700]. When the integration time is measured in units 
of the (Bondi) mass ms(0) on the initial slice, the last number is more than 600 000 m^(0). 
In figure ^| we plot the Bondi time u of an observer initially placed at the north pole of J 
for the conformal time interval [0, 1]. Most conspicuous is the rapid growth near i + , with 
the largest value u ~ 950. 

C. The decay of the Bondi mass 

Here we shortly sketch how we calculate the Bondi mass. Details are given in M. 
The Bondi mass of a cut S u of J is given by 

m B — — (V>2 + frfrj d 2 S, (6) 




0.0 



: * 












: + 












: x 



























































































































10 



Proper time of observer near i + 
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FIG. 8: Proper time along an observer world 
line near i + for a 50 3 (+), an 100 3 (x), and an 
150 3 (*) run. 
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FIG. 9: Bondi time of the Bondi observer at the 
north pole as a function of conformal time. 



where ~ denotes quantities with respect to a conformal metric g a b = o:~ 2 g a b in which J 
is expansion free. The quantity a is one of the spin coefficients in the Newman-Penrose 
formalism, its Bondi time derivative a is the news function. 

When calculating the world lines of Bondi observers, we keep track of the evolution of a for 
each observer and we parallelly transport with respect to g a b the null frame associated with 
the observer. 

Since the u = const cuts do in general not coincide with the t = const cuts, we have to locate 
the u = const cuts and interpolate onto them from the t = const cuts before we perform the 
integration. 
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Once the initial value of the Bondi mass mg(0) is determined, there is an alternative calcu- 
lation of the Bondi mass, using the mass loss formula: 

ms(M)=mB(0)- / / (&& d 2 S) du . (7) 
Jo Js u , v ' 

The latter method usually gives much smaller errors. Hence we use (^) to calculate the time 
evolution of the Bondi mass. The initial value is taken from the most accurate calculation, 
namely the 150 3 run. 

In our setup a is on the initial slice. Figure [lOj shows the integrand ip2 on the initial cut of 




FIG. 10: tp2 on the initial cut of J . Values range 
from -0.148 (black) to 0.115 (white). 



J . Obviously there is a strong high order moment. Since the maximal value of is about 
0.15, and the initial value of the Bondi mass is about 0.00104, the high order moments are 
by a factor of order 100 stronger than the monopole moment, which is the Bondi mass. 



Figure [XT] shows the Bondi mass as a function of the Bondi time u on a logarithmic scale 
(the abscissa is log (it + 0.1)). We observe a rapid decay in a first stage which lasts until 
u m 0.15. Then there is a slower decay lasting until u ~ 1.1. 

To see what happens after this stage we change the scale of the ordinate in figure where 
we plot the Bondi mass against a non-logarithmic time scale. 

Since the Bondi mass has already decayed by a factor of about 200, even our small numerical 
error becomes an issue. The result from the 50 3 (indicated by x) run significantly differs 
from the results for the 100 3 (+) run, whereas the later almost coincides with the results of 
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FIG. 11: Bondi mass as a function of Bondi 
time. 



3ondi mass at late times 
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FIG. 12: Bondi mass as a function of Bondi time 
at late times for 50 3 (x), 100 3 (+), and 150 3 (*) 
runs. 



the 150 3 (*) run. 

Analytically the Bondi mass of a spacetime with a regular i + must vanish at i + . It is not 
clear, whether the final value of the computed Bondi mass would eventually vanish, if we 
integrated even longer, or whether this offset is due to a numerical error in the initial value 
for the Bondi mass, which mainly depends on how accurate the provided initial data solve 
the constraints. 
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IV. CONCLUSION 

We have calculated the future time evolution of hyperboloidal gravitational wave data, which 
do not possess any continuous symmetry. We have seen that the conformal approach allows 
us to cover the entire physical future of these data with a finite grid and to determine 
the decay of curvature invariants over ranges unreachable by codes working in physical 
spacetime. 

Due to the use of higher order methods we can use fairly coarse grids. Since the grids used 
are coarse, we only need moderate amounts of computer resources, in particular our time 
evolution on a Origin2000 with R10000 processors requires 

• less than 15 minutes on 8 processors for a 50 3 run, where we need 87 time steps to 
cover i + , 

• less than 2 hours on 16 processors for a 100 3 run, where we need 173 time steps to 
cover i + , and 

• less than 6 hours on 27 processors for a 150 3 run, where we need 263 time steps to 
cover i + . 

Already in the 50 3 run the error is at most a few percent. With the 150 3 run we achieve an 
error of less than one part in thousand. 
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